arXiv:1509.04682v2 [math.OC] 7 Nov 2015 


Robust Sensitivity Analysis of the 
Optimal Value of Linear Programming 

Guanglin Xu* Samuel Bureh 

September 14, 2015 
Revised: November 4, 2015 


Abstract 

We propose a framework for sensitivity analysis of linear programs (LPs) in minimiza¬ 
tion form, allowing for simultaneous perturbations in the objective coefficients and 
right-hand sides, where the perturbations are modeled in a compact, convex uncer¬ 
tainty set. This framework unifies and extends multiple approaches for LP sensitiv¬ 
ity analysis in the literature and has close ties to worst-case linear optimization and 
two-stage adaptive optimization. We define the minimum (best-case) and maximum 
(worst-case) LP optimal values, p~ and p'^, over the uncertainty set, and we discuss 
issues of finiteness, attainability, and computational complexity. While p~ and are 
difficult to compute in general, we prove that they equal the optimal values of two 
separate, but related, copositive programs. We then develop tight, tractable conic 
relaxations to provide lower and upper bounds on p~ and p'^, respectively. We also 
develop techniques to assess the quality of the bounds, and we validate our approach 
computationally on several examples from—and inspired by—the literature. We find 
that the bounds on p~ and p'^ are very strong in practice and, in particular, are at 
least as strong as known results for specific cases from the literature. 

Keywords: Sensitivity analysis, minimax problem, nonconvex quadratic programming, 
semidefinite programming, copositive programming, uncertainty set. 


1 Introduction 

The standard-form linear program (LP) is 

min c^x 

s. t. Ax = b (1) 

a; > 0 
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where x G M"" is the variable and {A, b, c) G x x M"" are the problem parameters. 

In practice, {A, b, c) may not be known exactly or may be predicted to change within a 
certain region. In such cases, sensitivity analysis (SA) examines how perturbations in the 
parameters affect the optimal value and solution of Q. Ordinary SA considers the change of 
a single element in {A, b, c) and examines the corresponding effects on the optimal basis and 
tableau; see |8]. SA also extends to the addition of a new variable or constraint, although 
we do not consider such changes in this paper. 

Beyond ordinary SA, more sophisticated approaches that allow simultaneous changes in 
the coefficients c or right-hand sides b have been proposed by numerous researchers. Bradley 
et al. |1] discuss the 100-percent rule that requires specihcation of directions of increase or 
decrease from each cj and then guarantees that the same basis remains optimal as long as 
the sum of fractions, corresponding to the percent of maximum change in each direction 
derived from ordinary SA, is less than or equal to 1. Wendell [351 EEl EH develops the 
tolerance approach to hnd the so-called maximum tolerance percentage by which the objective 
coefficients can be simultaneously and independently perturbed within a priori bounds. The 
tolerance approach also handles perturbations in one row or column of the matrix coefficients 
[3T] or even more general perturbations in all elements of the matrix coefficients under certain 
assumptions [3H- Freund [15] investigates the sensitivity of an LP to simultaneous changes 
in matrix coefficients. In particular, he considers a linear program whose coefficient matrix 
depends linearly on a scalar parameter 9 and studies the effect of small perturbations on the 
the optimal objective value and solution; see also [^T] [221 EH]- Readers are referred to [31] 
for a survey of approaches for SA of problem ([^. 

An area closely related to SA is interval linear programming (ILP), which can be viewed as 
multi-parametric linear programming with independent interval domains for the parameters 
[IHlIHlEe]. Steuer [33] presents three algorithms for solving LPs in which the objective 
coefficients are specihed by intervals, and Gabrel et al. [I6] study LPs in which the right- 
hand sides vary within intervals and discuss the maximum and minimum optimal values. 
Mraz m considers a general situation in which the matrix coefficients and right-hand sides 
change within intervals and calculates upper and lower bounds for the associated optimal 
values. A comprehensive survey of ILP has been given by Hladik I2D1. 

To the best of our knowledge, in the context of LP, no authors have considered simulta¬ 
neous LP parameter changes in a general way, i.e., perturbations in the objective coefficients 
c, right-hand sides b, and constraint coefficients A within a general region (not just inter¬ 
vals). The obstacle for doing so is clear: general perturbations lead to nonconvex quadratic 
programs (QPs), which are NP-hard to solve (as discussed below). 

In this paper, we extend—and in many cases unify—the SA literature by employing 
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modern tools for nonconvex QPs. Specifically, we investigate SA for LPs in which (b, c) 
may change within a general compact, convex set U, called the uncertainty set. Our goal 
is to calculate—or bound—the corresponding minimum (best-case) and maximum (worst- 
case) optimal values. Since these values involve the solution of nonconvex QPs, we use 
standard techniques from copositive optimization to reformulate these problems into convex 
copositive programs (COPs), which provide a theoretical grounding upon which to develop 
tight, tractable convex relaxations. We suggest the use of semidefinite programming (SDP) 
relaxations, which also incorporate valid conic inequalities that exploit the structure of the 
uncertainty set. We refer the reader to [10] for a survey on copositive optimization and 
its connections to semidehnite programming. Relevant dehnitions and concepts will also be 


given in this paper; see Section 1.1 


Our approach is related to the recent work on worst-case linear optimization introduced 
by Peng and Zhu [29| in which: (i) only b is allowed to change within an ellipsoidal region; and 
(ii) only the worst-case LP value is considered. (In fact, one can see easily that, in the setup 
of [29] based on (i), the best-case LP value can be computed in polynomial time via second- 
order-cone programming, making it less interesting to study in their setup.) The authors 
argue that the worst-case value is NP-hard to compute and use a specialized nonlinear 
semidehnite program (SDP) to bound it from above. They also develop feasible solutions to 
bound the worst-case value from below and show through a series of empirical examples that 
the resulting gaps are usually quite small. Furthermore, they also demonstrate that their 
SDP-based relaxation is better than the so-called affine-rule approximation (see 0) and the 
Lasserre linear matrix ineguality relaxation (see pTl IT9] ). 

Our approach is more general than [29] because we allow both b and c to change, we 
consider more general uncertainty sets, and we study both the worst- and best-case values. 
In addition, instead of developing a specialized SDP approach, we make use of the machinery 
of copositive programming, which provides a theoretical grounding for the construction of 
tight, tractable conic relaxations using existing techniques. Nevertheless, we have been 
inspired by their approach in several ways. For example, their proof of NP-hardness also 
shows that our problem is NP-hard; we will borrow their idea of using primal solutions to 
estimate the quality of the relaxation bounds; and we test some of the same examples. 

We mention two additional connections of our approach with the literature. In [3], 
Bertsimas and Goyal consider a two-stage adaptive linear optimization problem under right- 
hand side uncertainty with a min-max objective. A simplihed version of this problem, in 
which the hrst-stage variables are non-existent, reduces to worst-case linear optimization; 
see the introduction of [3]. In fact, Bertsimas and Goyal use this fact to prove that their 
problem is NP-hard via the so-called max-min fractional set cover problem, which is a specihc 
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worst-case linear optimization problem studied by Feige et al. [13]. Our work is also related 
to the study of adjustable robust optimization |2lE9], which allows for two sets of decisions— 
one that must be made before the uncertain data is realized, and one after. In fact, our 
problem can viewed as a simplified case of adjustable robust optimization having no first- 
stage decisions. On the other hand, our paper is distinguished by its application to sensitivity 
analysis and its use of copositive and semidefinite optimization. 

We organize the paper as follows. In Section]^ we extend many of the existing approaches 
for SA by considering simultaneous, general changes in (6, c) and the corresponding effect 
on the LP optimal value. Precisely, we model general perturbations of (6, c) within a com¬ 
pact, convex set U —the uncertainty set, borrowing terminology from the robust-optimization 
literature—and define the corresponding minimum and maximum optimal values p~ and 
respectively. We call our approach robust sensitivity analysis, or RSA. Then, continuing in 
Section we formulate the calculation of p~ and p'^ as nonconvex bilinear QPs (or BQPs) 
and briefly discuss attainability and complexity issues. We also discuss how p~ and p'^ may 
be infinite and suggest alternative bounded variants, q~ and q'^, which have the property 
that, if p~ is already hnite, then q~ = p~ and similarly for g’*' and p^. Compared to re¬ 
lated approaches in the literature, our discussion of finiteness is unique. We then discuss the 
addition of redundant constraints to the formulations of q~ and g"*", which will strengthen 
later relaxations. Section then establishes COP reformulations of the nonconvex BQPs by 
directly applying existing reformulation techniques. Then, based on the COPs, we develop 
tractable SDP-based relaxations that incorporate the structure of the uncertainty set U, and 
we also discuss procedures for generating feasible solutions of the BQPs, which can also be 
used to verify the quality of the relaxation bounds. In Section]^ we validate our approach on 
several examples, which demonstrate that the relaxations provide effective approximations 
of g"*" and q~. In fact, we find that the relaxations admit no gap with g"*" and q~ for all tested 
examples. 

We mention some caveats about the paper. First, we focus only on how the optimal 
value is affected by uncertainty, not the optimal solution. We do so because we believe 
this will be a more feasible first endeavor; determining how general perturbations affect 
the optimal solution can certainly be a task for future research. Second, as mentioned 
above, we believe we are the hrst to consider these types of general perturbations, and thus 
the literature with which to compare is somewhat limited. However, we connect with the 
literature whenever possible, e.g., in special cases such as interval perturbations and worst- 
case linear optimization. Third, since we do not make any distributional assumptions about 
the uncertainty of the parameters, nor about their independence or dependence, we believe 
our approach aligns well with the general sprit of robust optimization. It is important to 
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note, however, that our interest is not robust optimization and is not directly comparable to 
robust optimization. For example, while in robust optimization one wishes to hnd a single 
optimal solution that works well for all realizations of the uncertain parameters, here we are 
only concerned with how the optimal value changes as the parameters change. Finally, we 
note the existence of other relaxations for nonconvex QPs including LP relaxations (see [3^ ) 
and Lasserre-type SDP relaxations. Generally speaking, LP-based relaxations are relatively 
weak (see ID); we do not consider them in this paper. In addition, SDP approaches can often 
be tailored to outperform the more general Lasserre approach as has been demonstrated in 
[29]. Our copostive- and SDP-based approach is similar; see for example the valid inequalities 
discussed in Section |3(2l 

1.1 Notation, terminology, and copositive optimization 

Let M” denote n-dimensional Euclidean space represented as column vectors, and let M” 
denote the nonnegative orthant in M"’. For a scalar p > 1, the p-norm of v G M” is dehned 
ll'^llp •= (X^r=i 6-g-! Il'^lli = Id|- We will drop the subscript for the 2-norm, 

e.g., ||n|| := ||ni| 2 - For v,w E M”, the inner product of v and w is dehned as v^w = YYi=i 
and the Hadamard product of v and w is dehned hj v o w := (vitci,..., G M"'. 

^mxn (denotes the set of real m x n matrices, and the trace inner product of two matrices 
A,Be is dehned A • B := trace(y4^i?). S"' denotes the space of n x n symmetric 

matrices, and for X G 5”, X ^ 0 denotes that X is positive semidehnite. In addition, 
diag(X) denotes the vector containing the diagonal entries of X. 

We also make several dehnitions related to copositive programming. The nxn copositive 
cone is dehned as 


COV{Wl) := {M eS^ : x^Mx > 0 V a: G }, 

and its dual cone, the completely positive cone, is 

CV{Wl) :={X eS^:X = e Rl}, 

where the summation over k is hnite but its cardinality is unspecihed. The term copositive 
programming refers to linear optimization over COV{W\_) or, via duality, linear optimization 
over CV{R^). A more general notion of copositive programming is based on the following 
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ideas. Let /C C M"- be a closed, convex cone, and define 

COVilC) := {M : x^Mx > 0 V x G /C}, 

CVilC) ■={XeS^ ■.X = e /C}. 

Then generalized copositive programming is linear optimization over COV{lC) and CV{X) 
and is also sometimes called set-semidefinite optimization m- In this paper, we work 
with generalized copositive programming, although we will use the shorter phrase copositive 
programming for convenience. 

2 Robust Sensitivity Analysis 

In this section, we introduce the concept of robust sensitivity analysis of the optimal value of 
the linear program Q. In particular, we define the best-case optimal value p~ and the worst- 
case optimal value p'^ over the uncertainty set U, which contains general perturbations in 
the objective coefficients c and the right-hand sides b. We then propose nonconvex bilinear 
QPs (BQPs) to compute p~ and p'^. Next, we clarify when p~ and could be infinite 
and propose finite, closely related alternatives and q~, which can also be formulated 
as nonconvex BQPs. Importantly, we prove that q~ equals p~ whenever p~ is finite; the 
analogous relationship is also proved for g"*" and p"*". 

2.1 The best- and worst-case optimal values 

In the Introduction, we have described h and c as parameters that could vary, a concept 
that we now formalize. Hereafter, (6, c) denotes the nominal, “best guess” parameter values, 
and we let (6, c) denote perturbations with respect to (6, c). In other words, the true data 
could be (h + h,c + c), and we think of h and c as varying. We also denote the uncertainty 
set containing all possible perturbations {b, c) as W C M™- x M”. Throughout this paper, we 
assume the following: 

Assumption 1. U is compact and convex, andU contains (0,0). 

Given (6, c) G U, we define the perturbed optimal value function at (6, c) as 

p{b, c) := min{(c -|- c)^x : Ax = b + b, x > 0}. (2) 

For example, p(0, 0) is the nominal optimal value of the nominal problem based on the nomi¬ 
nal parameters. The main idea of robust sensitivity analysis is then to compute the infimum 
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(best-case) and supremum (worst-case) of all optimal values p{b, c) over the uncertainty set 
U, i.e., to calculate 


p := inf{p(6, c) : (6, c) G W}, 
p'^ := sup{p(6, c) : (6, c) G W}. 

We illustrate p~ and with a small example. 

Example 1. Consider the nominal LP 


min Xi + X 2 
s. t. Xi + X 2 = 2 

Xi,X2 > 0 


and the uncertainty set 


U := 



bi G [-1,1] 

Cl G [-0.5, 0.5], C2 



( 3 ) 

( 4 ) 


( 5 ) 


Note that the perturbed data bi + bi and ci + ci remain positive, while C 2 + C 2 is constant. 
Thus, the minimum optimal value p~ occurs when bi and Ci are minimal, i.e., when bi = —1 
and Cl = —0.5. In this case, p~ = 0.5 at the solution {xi,X2) = (1,0). In a related manner, 
p'^ = 3 when bi = 1 and ci = 0.5 at the point {xi,X 2 ) = (0,3). Actually, any perturbation 
with Cl G [0,0.5] and bi = 1 realizes the worst-case value p^ = 3. Figure^ illustrates this 
example. 

We can obtain a direct formulation of p~ by simply collapsing the inner and outer mini¬ 
mizations of ([^ into a single nonconvex BQP: 

(c + cYx 

s. t. Ax = b-\-b, x>0 (6) 

(&, c) G U. 

The nonconvexity comes from the bilinear term c^x in the objective function. In the special 
case that (5, c) G U implies c = 0, i.e., when there is no perturbation in the objective 
coefficients, we have the following: 

Remark 1. If U is tractable and c = 0 for all {b,c) G U, then p~ can he computed in 
polynomial time as the optimal value of & with c = 0, which is a convex program. 
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(a) Illustration of the best-case optimal value 



(b) Illustration of the worst-case optimal value 


Figure 1: Illustration of Example Note that the dashed line in both (la) and (lb) 
corresponds to the feasible region of the nominal problem. 


A direct formulation for p~^ can, under a fairly weak assumption, be gotten via duality. 
Define the perturbed primal and dual feasible sets for any ( 6 , c) eU: 

P{b) := {x : Ax = b + b,x > 0}, 

D{c) := {{y,s) > 0 : A^y -|- s = c -|- c, s > 0}. 

For instance, P(0) and D(0) are the primal-dual feasible sets of the nominal problem. Next 
define the dual LP for (|^ as 

d{b, c) := max {(6 -|- b)^y : {y, s) G D{c)}. 

Considering the extended notion of strong duality, which handles the cases of infinite values, 
we have that d{b, c) = p{b, c) when at least one of P{b) and D{c) is nonempty. Hence, under 
the assumption that every ( 6 , c) yields P{b) 7 ^ 0 or D{c) 7 ^ 0, a direct formulation for 
can be constructed by replacing p( 6 , c) in ([^ with d{b^ c) and then collapsing the subsequent 
inner and outer maximizations into the single nonconvex BQP 

s^Pb,c,y,s {b + bfy 

s. t. A'^y + s = c + c, s >0 (7) 

(6, c) e U. 


Here again, the nonconvexities arise due to the bilinear term b^y in the objective. If ( 6 , c) eU 
implies 6 = 0 , then p~^ can be catenated in polynomial time: 







































Remark 2. If U is tractable and 6 = 0 for all {b, c) G U, then can be computed in 
polynomial time as the optimal value of & with 6 = 0, which is a convex program. 

We summarize the above discussion in the following proposition: 

Proposition 1. The best-case value p~ equals the optimal value of 0. Moreover, ifP{b) ^ 0 
or D{c) 7 ^ 0 for all ( 6 , c) G U, then the worst-case value p'^ equals the optimal value of 

We view the condition in Proposition [I]— that at least one of -P( 6 ) and D{c) is nonempty 
for each ( 6 , c) G U —to be rather mild. Said differently, the case that P( 6 ) = D{c) = 0 for 
some ( 6 , c) eU appears somewhat pathological. For practical purposes, we hence consider 
Q to be a valid formulation of p'^. Actually, in the next subsection, we will further restrict 
our attention to those ( 6 , c) eU for which both P{b) and P{c) are nonempty. In such cases, 
each p{b, c) is guaranteed to be hnite, which—as we will show—carefully handles the cases 
when and p~ are inhnite. 

Indeed, the worst-case value p~^ could equal -|-oo due to some perturbed P( 6 ) being empty 
as shown in the following example: 

Example 2. In Example\^ change the uncertainty set to 

6 i G [—3,1] 

Cl G [-0.5,0.5], C2 = 0 

Then p{b, c) = -|-cxd whenever 6 i G [—3, —2) since then the primal feasible set P{h) is empty. 
Then p'^ = -|-cxd overall. However, limiting hi to [—2,1] yields a worst-case value of 3 as 
discussed in Example\^ 

Similarly, p~ might equal —oo due to some perturbed LP having unbounded objective value, 
implying infeasibility of the corresponding dual feasible set D{c). 

2.2 Attainability and complexity 

In this brief subsection, we mention results pertaining to the attainability of p~ and p'^ and 
the computational complexity of computing them. 

By an existing result concerning the attainability of the optimal value of nonconvex 
BQPs, we have that p~ and p"*" are attainable when U has a relatively simple structure: 

Proposition 2 (theorem 2 of [25]). Suppose U is representable by a finite number of linear 
constraints and at most one convex quadratic constraint. Then, if the optimal value of 0 
is finite, it is attained. A similar statement holds for 0 . 
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In particular, attainability holds when lA is polyhedral or second-order-cone representable 
with at most one second-order cone. Moreover, the bilinear nature of ([^ implies that, if the 
optimal value is attained, then there exists an optimal solution {x*,b*,c*) with {b*,c*) an 
extreme point of U. The same holds for Q if its optimal value is attained. 

As discussed in the Introduction, the worst-case value has been studied by Peng and 
Zhu [29] for the special case when c = 0 and b is contained in an ellipsoid. The authors 
demonstrate (see their proposition 1.1) that calculating in this case is NP-hard. By the 
symmetry of duality, it thus also holds that p~ is NP-hard to compute in general. 


2.3 Finite variants of p and p'^ 

We now discuss closely related variants of and p~ that are guaranteed to be hnite and 
to equal p'^ and p~, respectively, when those values are themselves hnite. We require the 
following feasibility and boundedness assumption: 

Assumption 2. Both feasible sets P{0) and D{0) are nonempty, and one is bounded. 

By standard theory, P(0) and hl(0) cannot both be nonempty and bounded. Also dehne 


U-.= {{b,c)eU-.P{b)^%,D{c)^%}. 


Note that (0, 0) G W due to Assumption]^ In fact, U can be captured with linear constraints 
that enforce primal-dual feasibility and hence is a compact, convex subset of U\ 


U 


{b, c) eU 


Ax = 
A^y -\- s 


b + b,x >t) 

= c -f c, s > 0 


Analogous to p'^ and p , dehne 


:= sup{p(6, c) : (6, c) eU] (8) 

q~ ;= mi{p{b, c) : (b, c) eU}. (9) 

The following proposition establishes the hniteness of q'^ and q~\ 

Propositions. Under Assumptions\^and^ both q^ and q~ are finite. 

Proof. We prove the contrapositive for q~. (The argument for q~^ is similar.) Suppose 
q~ = —oo. Then there exists a sequence {{b^, c^)} PU with hnite optimal values p{b^, c^) —)■ 
—oo. By strong duality, there exists a primal-dual solution sequence {{x^,y^,s^)} with 
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(c + c)'^= {h + hY'y^ —)■ —oo. Since U is bounded, it follows that ||a:^|| —)■ oo and 
\\y^\\ -)■ cx). 

Consider the sequence with := {x’^,b^)/\\x'^\\. We have Az'^ = 6/||a;^|| + 

d^, z^ > 0, and ||z^|| = 1 for all k. Morover, 6/||a;^|| + —?■ 0. Hence, there exists a 
subsequence converging to (^,0) such that Az = 0, ^ > 0, and ||f|| = 1. This proves that 
the recession cone of -P(O) is nontrivial, and hence -P(O) is unbounded. In a similar manner, 
D{0) is unbounded, which means Assumption does not hold. □ 

Note that the proof of Proposition only assumes that W, and hence W, is bounded, which 
does not use the full power of Assumption 

Similar to , a direct formulation of q~ can be constructed by employing the primal-dual 
formulation of U and by collapsing the inner and outer minimizations of (|^ into a single 
nonconvex BQP: 


q~ = mib,c,x,y,s {c + c)'^x 

s. t. Ax = b + b, a: > 0 

A^y -\- s = c + s>0 

(6, c) G U. 


( 10 ) 


Likewise for after replacing p{b, c) in (|^ by d(6, c), we can collapse the inner and outer 
maximizations into a single nonconvex BQP: 


g+ = {b + bfy 

s. t. Ax = b + b, X >0 

A^y + s = c + c, s>0 
(6, c) G U. 


( 11 ) 


The following proposition establishes q^ = when is hnite and, similarly, q = p 
when p~ is hnite. 

Proposition 4. If p'^ is finite, then q'^ = P'^, and if p~ is finite, then q~ = p~ ■ 


Proof. We prove the second statement only since the hrst is similar. Comparing the formu¬ 
lation (|^ for p~ and the formulation (10) for q~, it is clear that p~ < q~ ■ In addition, let 
(6, c, x) be any feasible solution of 


Because p~ is hnite, p{b, c) is hnite. Then the cor¬ 
responding dual problem is feasible, which implies that we can extend (6, c, x) to a solution 
(b,c,x,y, s) of (10) with the same objective value. Hence, p~ > q~. □ 

In the remaining sections of the paper, we will focus on the hnite variants q~ and q^ given 
by the nonconvex QPs (10) and which optimize the optimal value function p{b, c) = 
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d{b, c) based on enforcing primal-dual feasibility. It is clear that we may also enforce the 
complementary slackness condition x o s = 0 without changing these problems. Although 
it might seem counterintuitive to add the redundant, nonconvex constraint x o s = 0 to an 
already difficult problem, in Section we will propose convex relaxations to approximate 
q~ and q~^, in which case—as we will demonstrate—the relaxed versions of the redundant 
constraint can strengthen the relaxations. 


3 Copositve Formulations and Relaxations 


In this section, we use copositive optimization techniques to reformulate the RSA problems 
(10) and (0 into convex programs. We further relax the copositive programs into conic, 
SDP-based problems, which are computationally tractable. 


3.1 Copositive formulations 


In order to formulate (10) and ( [II] ) as COPs, we apply a result of [6]; see also [5] [O] [12]. 
Consider the general nonconvex QP 


inf z'^Wz + 2w'^z 
s. t. Ez = /, z & K, 


( 12 ) 


where /C is a closed, convex cone. Its copositive reformulation is 

inf W • Z + 2 z 
s. t. Ez = /, diag(EZE^) = / o / 
1 z'^" 


z Z 


G CP(M+ X /C), 


as established by the following lemma: 


(13) 


Lemma 1 (corollary 8.3 in [6]). Problem (12) is equivalent to (13), i.e.: (i) both share the 


same optimal value; (ii) if {z*, Z*) is optimal for (13), then z* is in the convex hull of optimal 


solutions for (12). 


The following theorem establishes that problems (10) and ( [II| can be reformulated as 
copositive programs according to Lemma [Tj The proof is based on describing how the two 


problems £t the form (12). 
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Theorem 1. Problems (10) and (11) to compute q and g"* 
grams of the form (13), where 


are solvable as copositive pro- 


K. ;= hom(W) x 


X X 


and 


X X 


:t>0, (6,c)/teW}U {(0,0,0)} 


hom(W) := {{t,h,c) G 
is the homogenization ofU. 

Proof. We prove the result for just problem ( 10 ) since the argument for problem (|IT| 


IS 


similar. First, we identify z G /C in (12) with {t,b,c,x,y,s) G hom(W) x 


(10). In addition, in the constraints, we identify Ez = f with the equations Ax = tb + b, 


A'^y + s = tc + c, and t = 1. Note that the right-hand-side vector / is all zeros except for 
a single entry corresponding to the constraint t = 1. Moreover, in the objective, z'^Wz is 
identified with the bilinear term c^x, and 2w'^z is identified with the linear term c^x. With 
this setup, it is clear that (10) is an instance of (12) and hence Lemma applies to complete 
the proof. □ 


3.2 SDP-based conic relaxations 


As discussed above, the copositive formulations of ( 10 ) and 0 as represented by ( [I^ are 
convex yet generally intractable. Thus, we propose SDP-based conic relaxations that are 
polynomial-time solvable and hopefully quite tight in practice. In Section below, we will 
investigate their tightness computationally. 


We propose relaxations that are formed from (13) by relaxing the cone constraint 


M := I ^ ^ 1 G CP(M+ X /C). 
X Z 


As is well known—and direct from the dehnitions—cones of the form CV{-) are contained 
in the positive semidefinite cone. Hence, we will enforce M ^ 0. It is also true that 
M G C'P(M+ X /C) implies z & IC, although M P 0 does not necessarily imply this. So, in 
our relaxations, we will also enforce z E K,. Including z E K, improves the relaxation and 
also helps in the calculation of bounds in Section |3.3| 

Next, suppose that the description of M+ x /C contains at least two linear constraints, 
afz < hi and a^z < b 2 . By multiplying bi — afz and 62 — ^ 2 ^, we obtain a valid, yet 
redundant, quadratic constraint 6162 — bia^z — b 2 afz + afzz'^a 2 > 0 for CP(]R+ x /C). This 
quadratic inequality can in turn be linearized in terms of M as hih 2 — hia"^z — h 2 afz-{-afZa 2 > 
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0, which is valid for CP(]R+ x /C). We add this linear inequality to our relaxation; it is called 
an RLT constraint |32]. In fact, we add all such RLT constraints arising from all pairs of 
linear constraints present in the description of M+ x K,. 

When the description of ]R_|_ x /C contains at least one linear constraint ajz < bi and 
one second-order-cone constraint \\d 2 — C^zW < 62 — ci^z, where ^2 is a vector and C 2 is a 
matrix, we will add a so-called SOC-RLT constraint to our relaxation [7j. The constraint is 
derived by multiplying the two constraints to obtain the valid quadratic second-order-cone 
constraint 

||( 6 i - alz){d 2 - C 2 z)\\ < {hi - alz){h 2 - a^z). 

After linearization by M, we have the second-order-cone constraint 

\\b 1 d 2 — d 2 aj'z — biC 2 z + C 2 Zai\\ < 6162 — bia^z — b 2 a\z -\- a\Za 2 - 


Finally, recall the redundant complementarity constraint xos = 0 described at the end of 
Section 2.3, which is valid for both (10) and ( pT| . Decomposing it as XiSi = 0 for i = 1,..., n, 
we may translate these n constraints to ( |T^ as z^HiZ = 0 for appropriatly dehned matrices 
matrices Hi. Then they may be linearized and added to our relaxation as • Z = 0. 

To summarize, let RLT denote the set of {z, Z) satisfying all the derived RLT constraints, 
and similarly, dehne SOCRLT as the set of {z, Z) satisfying all the derived SOC-RLT con¬ 
straints. Then the SDP-based conic relaxation for (13) that we propose to solve is 


inf 

s. t. 


W • Z + 2w'^z 


Ez = /, diag{EZE^) = fof 
Hi • Z = 0 Vi = l,...,n 
{z, Z) e RLT n SOCRLT 




^ 0, Z e /C. 


(14) 


It is worth mentioning that, in many cases, the RLT and SOC-RLT constraints will already 


imply z E }C, but in such cases, we nevertheless write the constraint in (14) for emphasis; 
see also Section [3731 below. 


When translated to the problem (10) for calculating q , the relaxation (14) gives rise to 
a lower bound < q~. Similarly, when applied to (11), we get an upper bound g^p > g+. 
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3.3 Bounds from feasible solutions 


In this section, we discuss two methods to approximate q~ from above and from below, 
i.e., to bound q~ and q~^ using feasible solutions of (10) and ( [TT| ), respectively. 

The first method, which has been inspired by [29], utilizes the optimal solution of the 


SDP relaxation (14). Let us discuss how to obtain such a bound for (10), as the discussion 


for ([II| is similar. We first observe that any feasible solution {z, Z) of (14) satisfies Ez = / 


and z ^ 1C, i.e., 2 ; satisfies all of the constraints of (12). Since (12) is equivalent to (10) 


under the translation discussed in the proof of Theorem z gives rise to a feasible solution 
{x, y, s, b, c) of (10). From this feasible solution, we can calculate (c + c)^x > q~. In practice, 

z- 


we will start from the optimal solution {z 
following remark. 


of (14). We summarize this approach in the 


Remark 3. Suppose that {z , Z ) is an optimal solution of the relaxation {I 4 ) corresponding 


to and let {x ,y ,s ,b ,c ) be the translation of z to a feasible point of ([i^. Then, 
r~ := (c + c~Yx~ > q~. Similarly, we define r~^ := {b + < q'^ based on an optimal 

solution {z~^,Z~^) of { 14 ) corresponding to 0. 


Our second method for bounding q~ and g"*" using feasible solutions is a sampling proce¬ 
dure detailed in Algorithm The main idea is to generate randomly a point (6, c) eU and 
then to calculate p{b,c), which serves as an upper bound of p~ and a lower bound of 
i.e., p~ < p{b,c) < p'^. Multiple points {b^,c^) and values p^ := p{b^,c^) are generated and 
the best bounds p~ < v~ := min^jp^} and max^jp^} =: v~^ < p'^ are saved. In fact, by the 


bilinearity of (10) and (11), we we may restrict attention to the extreme points (&, c) of U 


without reducing the quality of the resultant bounds; see also the discussion in Section |2.2 
Hence, Algorithm generates—with high probability—a random extreme point of U by op¬ 
timizing a random linear objective over U, and we generate the random linear objective as a 
vector uniform on the sphere, which is implemented by a well-known, quick procedure. Note 
that, even though the random objective is generated according to a specific distribution, we 
cannot predict the resulting distribution over the extreme points of U. 

As all four of the bounds r~,r~^,v~, and n"*" are constructed from feasible solutions, we 


can further improve them heuristically by exploiting the bilinear objective functions in (10) 


and ( pT| . In particular, we employ the standard local improvement heuristic for programs 
with a bilinear objective and convex constraints (e.g., see [23]). Suppose, for example, that 
we have a feasible point {x~ ,y~, s~ ,b~ ,c~) for problem (10) as discussed in Remark]^ To 


attempt to improve the solution, we fix the variable c in (10) at the value c , and we solve the 


resulting convex problem for a new, hopefully better point (x^, y^, s^,E, c^), where = c . 
Then, we fix x to x^, resolve, and get a new point (x^, s^, 6^, c^), where x^ = x^. This 
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Algorithm 1 Sampling procedure to bound q from above and from below 

Inputs: Instance with uncertainty set U and restricted uncertainty set U. Number of 
random trials T. 

Outputs: Bounds v~ := min^jp^} > p~ and := max^jp*^} < p"*". 

for k = 1,..., T do 

Generate (/, g) G x uniformly on the unit sphere. 

Calculate ( 6 ^, c^) G Argmin {/'^6 + g'^c : (6, c) G U}. 

Set p^ := c^). 

end for 


alternating process is repeated until there is no further improvement in the objective of ( 10 ), 
and the hnal objective is our bound r~. 

In Section 1^ below, we use the bounds r“, r+, and to verify the quality of our 
bounds and g^p. Our tests indicate that neither bound, r~ nor dominates the 
other—and similarly for the bounds and . Hence, we will actually report the better of 
each pair: min{r“,u“} and max{r’'', u’*'}. Also, for the calculations of v~ and u’*', we always 
take T = 10, 000 in Algorithmic 


4 Computational Experiments 

In this section, we validate our approach by testing it on six examples from the literature as 
well as an example of our own making. The hrst three examples in Section [4T] correspond to 
classical sensitivity analysis approaches for LP; the fourth example in Sectioncorresponds 
to an interval LP in inventory management; the hfth example in Section [4^ corresponds to 
a systemic-risk calculation in hnancial systems; and the last example in Section |4.4| is a 


transportation network flow problem. We implement our tests in Python (version 2.7.6) 
with Mosek (version 7.1.0.33) as our convex-optimization solver. All of Mosek’s settings are 
set at their defaults, and computations are conducted on a Macintosh OS X Yosemite system 
with a quad-core 3.20GHz Intel Gore i5 GPU and 8 GB RAM. 


4.1 Examples from classical sensitivity analysis 

Gonsider the following nominal problem from 


min — 12 xi — 18 x 2 — 18 x 3 — 40x4 
s. t. 4xi -|- 9 X 2 + 7 x 3 + 10 x 4 + 

X2 + 3X3 + 40X4 + Xq 

Xi,. . . ,X6 > 0. 


6000 

4000 
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The optimal basis is -B = {1,4} with optimal solution 1(4000,0,0,200,0,0) and optimal 
value p(0, 0) = —18667. According to standard, “textbook” sensitivity analysis, the optimal 
basis persists when the coefficient of xi lies in the interval [—16, —10] and other parameters 
remain the same. Along this interval, one can easily compute the best-case value —24000 
and worst-case value —16000, and we attempt to reproduce this analysis with our approach. 
So let us choose the uncertainty set 


r 6i = 62 = 0 1 

W = i (6, c) e X : Cl e [-4,2] > , 

[ C 2 = ■ ■ ■ = C6 = 0 J 

which corresponds precisely to the above allowable decrease and increase on the coefficient of 
xi- Note that Assumptions[T]and[^are satished. We thus know from above that = —24000 
and = —16000. Since 6 = 0 in W, Remark implies that q~^ is easy to calculate. So 
we apply our approach, i.e., solving the SDP-based relaxation, to approximate q~. The 
relaxation value is = —24000, which recovers q~ exactly. The CPU time for computing 
g^p is 0.10 seconds. 

Our second example is also based on the same nominal problem from [M], but we consider 
the 100%-rule. Again, we know that the optimal basis B = (1,4} persists when the coefficient 
of Xi lies in the interval [—16, —10] (and all other parameters remain the same) or separately 
when the coefficient of X 2 lies in the interval [—134/3, -|-oo] (and all other parameters remain 
the same). In accordance with the 100%-rule, we choose to decrease the coefficient of Xi, and 
thus its allowed interval is [—16, —12] of width 4. We also choose to decrease the coefficient of 
X 2 , and thus its allowed interval is [—134/3, —18] of width 80/3. The 100%-rule ensures that 
the optimal basis persists as long as the sum of fractions, corresponding to the percent of 
maximum changes in the coefficients of xi and X2, is less than or equal to 1. In other words, 
suppose that ci and C2 are the perturbed values of the coefficients of xi and 0:2, respectively, 
and that all other coefficients stay the same. Then the nominal optimal basis persists for 
(ci,C 2 ) in the following simplex: 


r Cl e [-16,-12] } 

<(ci, 52): C2 e [-134/3,-18] 1. 

I —12—Cl —18—C 2 ^ 1 I 

4 80/3 - 

By evaluating the three extreme points (-12,-18), (-16,-18) and (-12,-134/3) of this 
set with respect to the nominal optimal solution, one can calculate the best-case optimal 
value as q~ = —24000 and the worst-case optimal value as q'^ = —18667. We again apply 
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our approach in an attempt to recover empirically the 100%-rule. Specihcally, let 


U = 


6i = 62 = 0 

/, N Cl e [-4,0], C2 e [-y,o] 

4 80/3 — 

C3 = ■ ■ ■ = C6 = 0 


Note that Assumptions and are satished. Due to 6 = 0 and Remark we focus our 
attention on q~. Calculating the SDP-based relaxation value, we see that = —24000, 
which recovers q~ precisely. The CPU time is 0.15 seconds 

Our third example illustrates the tolerance approach, and we continue to use the same 
nominal problem from [M]. As mentioned in the Introduction, the tolerance approach consid¬ 
ers simultaneous and independent perturbations in the objective coefficients by calculating 
a maximum tolerance percentage such that, as long as selected coefficients are accurate to 
within that percentage of their nominal values, the nominal optimal basis persists; see [H7] . 
Let us consider perturbations in the coefficients of xi and X 2 with respect to the nominal 
problem. Applying the tolerance approach of EH. the maximum tolerance percentage is 1/6 
in this case. That is, as long as the two coefficient values vary within —12±12/6 = [—14, —10] 
and —18 ± 18/6 = [—21, —15], respectively, then the nominal optimal basis B = {1,4} per¬ 
sists. By testing the four extreme points of the box of changes [—14, —10] x [—21, —15] with 
respect to the optimal nominal solution, one can calculate the best-case optimal value as 
q~ = —21333 and the worst-case optimal value as q'^ = —16000. To test our approach in 
this setting, we set 


61 = 62 = 0, C3 = ■ ■ ■ = C6 = 0 I 
Cl G [— 2 , 2 ],C 2 G [—3,3] J 

and, as in the previous two examples, we focus on q~. Assumptions [I] and are again 
satished, and we calculate the lower bound g^p = —21333, which recovers q~ precisely. The 
CPU time for computing g^p is 0.13 seconds. 



4.2 An example from interval linear programming 

We consider an optimization problem that is typical in inventory management, and this 
particular example originates from [121 . Suppose one must decide the quantity to be ordered 
during each period of a hnite, discrete horizon consisting of T periods. The goal is to 
satisfy exogenous demands dk for each period fc, while simultaneously minimizing the total 
of purchasing, holding, and shortage costs. Introduce the following variables for each period 
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k: 


Sk = stock available at the end of period fc; 

Xk = quantity ordered at the beginning of period k. 


Items ordered at the beginning of period k are delivered in time to satisfy demand during 
the same period. Any excess demand is backlogged. Hence, each Xk is nonnegative, each Sk 
is free, and 

Sk—l d" Xk Sk dk- 


The order quantities Xk are further subject to uniform upper and lower bounds, u and /, and 
every stock level Sk is bounded above by U. At time k, the purchase cost is denoted as Ck, 
the holding cost is denoted as hk, and the shortage cost is denoted Qk- Then, the problem 
can be formulated as the following linear programming problem (assuming that the initial 
inventory is 0): 

+ Vk) 

s. t. So = 0 


Sfc— 1 T Xk Sk dk 

Vk ^ hkSk 
Uk ^ 9kSk 
I < Xk < U 
Sk <U 
Xk 5 Vk ^ 0 


fc = 1, . 
k = 1,. 
k = 1,. 
fc = 1,. 
k = 1,. 
fc = 1,. 


,T. 


(15) 


As in [16], consider an instance of (15) in which T = 4, u = 1500,/ = 1000, U = 600, 
and all costs are as in Table Moreover, suppose the demands dk are each uncertain 
and may be estimated by the intervals di G [700, 900], (i 2 ^ [1300,1600], da G [900,1100], 
and di G [500,700]. From [T6|, the worst-case optimal value over this uncertainty set is 
q~^ = 25600. For our approach, it is easy to verify that Assumptions [I] and are satisfied, 
and solving our SDP-based conic relaxation with an uncertainty set corresponding to the 
intervals on dk, we recover exactly, i.e., we have = 25600. The CPU time for 
computing our SDP optimal value is 1,542 seconds. 

Since the uncertainties only involve the right-hand sides. Remark implies that the 
best-case value q~ can be calculated in polynomial-time by solving an LP that directly 
incorporates the uncertainty. 
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Period {k) Purchasing cost (cfc) Holding cost {hk) Shortage cost {gk) 

1 7 2 3 

2 1 14 

3 10 1 3 

4 6 13 

Table 1: Costs for each period of an instance of the inventory management problem 

4.3 Worst-case linear optimization 

We next consider an example for calculating systemic risk in hnancial systems, which is an 
application of worst-case linear optimization (WCLO) presented in [2U] . 

For an interbank market, systemic risk is used to evaluate the potential loss of the whole 
market as a response to the decisions made by the individual banks [ISl EH]. Specihcally, 
let us consider a market consisting of n banks. We use an n x n matrix L to denote the 
liability relationship between any two banks in the market. For instance, the element 
represents the liability of bank i to bank j. In the market, banks can also receive exogenous 
operating cash flows to compensate their potential shortfalls on incoming cash flows. We use 
bi to denote the exogenous operating cash flow received by bank i. Given the vector b, we 
calculate the systemic loss l{b) of the market, which measures the amount of overall failed 
liabilities |29j : 


l(b) = min 

X 

s-t. {YTj=iLij)xi-YTj=iLjiXj <bi Vi = l,...,n 
Xi <1 V i = 1,..., n. 

Here the decision variables Xi represent the ratio of the total dollar payment by bank i to 
the total obligation of bank i. These ratios are naturally less than or equal to 1 < 1) as 

the banks do not pay more than their obligations. In contrast, 1 — Xi denotes the ratio of 
bank i failing to fulhll its obligations. Furthermore, we have a less-than-or-equal-to sign in 
the hrst constraint as the system allows limited liability (see [I3]). Finally, the objective is 
to minimize the total failure ratio of the system. 

In practice, however, there exist uncertainties in the exogenous operating cash flows. 
Allowing for uncertainties, the worst-case systemic risk problem [29] is given as 

max min (1 — xd 

feev a: 

s-t- {YJj=iLij)xi-YJj=iLjiXj <bi + Qi.b Vi = l,...,u 
Xi <1 V z = 1,..., n. 
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Instance 

max{r^, n’*'} 

gap+ 

?sdp 


1 

0.521 

0.0% 

0.521 

0.66 

2 

1.429 

0.0% 

1.429 

0.82 

3 

1.577 

0.0% 

1.577 

0.83 

4 

0.344 

0.0% 

0.344 

0.68 

5 

0.724 

0.0% 

0.724 

0.73 

6 

0.690 

0.0% 

0.690 

0.80 

7 

1.100 

0.0% 

1.100 

0.73 

8 

0.458 

0.0% 

0.458 

0.78 

9 

0.427 

0.0% 

0.427 

0.74 

10 

0.300 

0.0% 

0.300 

0.78 


Table 2: Results for the systemic-risk example 


where V := {6 G : ||6|| < 1} denotes the uncertainty set, Q G for some m < n 

corresponds to an affine scaling of V, and Qi. denotes the i-th row of Q. After converting 
the nominal LP to our standard form, we can easily put the systemic risk problem into our 
framework by dehning W := {(6, c) : 6 G V, c = 0} and slightly changing our lA to reflect the 
dependence on the affine transformation as represented by the matrix Q. 

Similar to table 4 in [29], we randomly generate 10 instances of size m x n = 3 x 5. In 
accordance with Remark [T| which states that q~ is easy to calculate in this case, we focus 
our attention on the worst-case value . It is straightforward to verify Assumptions [T] and 
In Table we list our 10 upper bound^ (one for each of the 10 instances) in the column 
titled and we report the computation time (in seconds) for all 10 instances under the 
column marked f^p. To evaluate the quality of g^p, we also calculate max{r+,n+} for each 
instance and the associated relative gap: 


gap+ 


Q'iir, ~ max{r+, 

- / + +11 n X 100^°- 

max{| max{r+,j|, 1| 


The computation times for computing r~ and are trivial, while the average computation 
time for computing v~ and n"'' is about 77 seconds. 

From the results in Table we see clearly that our approach recovers for all 10 
instances, which also matches the quality of results from [29] , 

^Mosek encountered numerical problems on some—but not all—of the generated system-risk instances. 
In Tablej^ we show 10 instances on which Mosek had no numerical issues. From private communication with 
the Mosek developers, it appears that the upcoming version of Mosek (version 8) will have fewer numerical 
issues on these instances. 
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Customer 




Supplier 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

1 





2 


3 


2 


2 

3 

3 

4 





1 

4 


3 


4 

5 

3 







4 

1 

2 


1 

3 

1 



8 

2 

5 

4 



3 



2 

1 

2 

1 


Table 3: The unit transportation costs associated with the arcs of the network. 

4.4 A network flow problem 

Next we consider a transportation network flow problem from |38], which has mi = 5 
suppliers/origins and m 2 = 10 customers/destinations for a total of m = 15 facilities. The 
network is bipartite and consists of n = 24 arcs connecting suppliers and customers; see 
Figure]^ Also shown in Figure]^ are the (estimated) supply and demand numbers (b) for each 
supplier and customer. In addition, the (estimated) unit transportation costs (c) associated 
with the arcs of the network are given in Table Suppose at the early stages of planning, 
the supply and demand units and the unit transportation costs are uncertain. Thus, the 
manager would like to quantify the resulting uncertainty in the optimal transportation cost. 



Figure 2: The transportation network of the 5 suppliers and 10 customers. 

We consider three cases for the uncertainty set, each of which is also parameterized by a 
scalar 7 G (0,1). In the hrst case (“POLY”), we consider the polyhedral uncertainty set 

^ 1 ( 7 ) = {(^c) : ||6||i < 7 || 6 ||i, ||c||i < 7 ||c||i }; 
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in the second case (“SOC”), we consider the second-order-cone uncertainty set 


^ 2 ( 7 ) := {{b,c) : || 6 || < 7 || 6 ||, ||c|| < 7 ||c|| }; 


and in the third case (“MIX”), we consider a mixture of the hrst two cases: 


^ 3 ( 7 ) := {{b,c) : ||6||i < 7 || 6 ||i, ||c|| < 7 ||c|| }. 


For each, 7 controls the perturbation magnitude in b and c relative to b and c, respectively. 
In particular, we will consider three choices of 7 : 0.01, 0.03, and 0.05. For example, 7 = 0.03 
roughly means that b can vary up to 3% of the magnitude of b. In total, we have three cases 
with three choices for 7 resulting in nine overall experiments. 

Assumptions and are satished in this example, and so we apply our approach to 
bound q~ and see Table Our 18 bounds (lower and upper bounds for each of the nine 
experiments) are listed in the two columns titled q~^p and respectively. We also report 
the computation times (in seconds) for all 18 instances under the two columns marked 
and t^p. We also compute r~, v~, r+, and v~^ and dehne the relative gaps 


gap = 
gap+ = 


min{r ,v } - 
max{| min{r“, n“}|, 1 } 

gs'dp -J^ax{r+,n+} 

max{| max{r+, n+}|, 1 } 


X 100%, 
X 100%. 


Again, the computation times for r~ and r~^ are trivial. The average computation time for 
computing v~ and is about 216 seconds. 

Table shows that our relaxations capture q~ and q~ in all cases. As ours is the hrst 
approach to study general perturbations in the literature, we are aware of no existing methods 
for this problem with which to compare our results. 


4.5 The effectiveness of the redundant constraint x o s = 0 


Finally, we investigate the effectiveness of the redundant complementarity constraint in (10) 
and (0 by also solving relaxations without the the linearized version of the constraint. As 
it turns out, in all calculations of q~^p, dropping the linearized complementarity constraint 
does not change the relaxation value. However, in all calculations of g^p, dropping it has a 
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Case 

7 

^sdp 

gap- 

min{r , v } 

maxjr’'', 

gap+ 


^sdp(s) 



0.01 

2638.4 

0.0% 

2638.4 

3088.8 

0.0% 

3088.8 

3551 

4259 

POLY 

0.03 

2139.6 

0.0% 

2139.6 

3437.4 

0.0% 

3437.4 

5106 

4433 


0.05 

1640.9 

0.0% 

1640.9 

3769.8 

0.0% 

3769.8 

5330 

4244 


0.01 

2745.6 

0.0% 

2745.6 

2981.6 

0.0% 

2981.6 

190 

122 

SOC 

0.03 

2498.9 

0.2% 

2504.3 

3212.1 

0.0% 

3212.1 

174 

112 


0.05 

2257.6 

0.2% 

2263.0 

3442.7 

0.0% 

3442.7 

150 

124 


0.01 

2724.1 

0.0% 

2724.1 

3008.4 

0.0% 

3008.4 

997 

838 

MIX 

0.03 

2429.2 

0.0% 

2429.2 

3281.9 

0.0% 

3281.9 

908 

731 


0.05 

2134.3 

0.0% 

2134.3 

3560.7 

0.0% 

3560.7 

1001 

851 


Table 4: Results for the transportation network problem 


Example_value without constraint 


Sec. 

4.1 

(#1) 

-16000 

100% 

0 

Sec. 

4.1 

(#2) 

-18667 

100% 

0 

Sec. 

4.1 

(#3) 

-16000 

100% 

0 

Sec. 

4.2 


25600 

1671% 

453298 

Sec. 

4.4 

(POLY 0.01) 

3088.8 

5.7% 

3265.8 

Sec. 

4.4 

(POLY 0.03) 

3437.4 

2.7% 

3528.5 

Sec. 

4.4 

(POLY 0.05) 

3769.8 

2.3% 

3855.6 

Sec. 

4.4 

(SOC 0.01) 

2981.6 

87.6% 

5593.1 

Sec. 

4.4 

(SOC 0.03) 

3212.1 

84.3% 

5920.2 

Sec. 

4.4 

(SOC 0.05) 

3442.7 

81.6% 

6252.7 

Sec. 

4.4 

(MIX 0.01) 

3008.4 

10.3% 

3319.4 

Sec. 

4.4 

(MIX 0.03) 

3281.9 

5.2% 

3453.5 

Sec. 

4.4 

(MIX 0.05) 

3560.7 

3.6% 

3689.4 

Table 5: 

Effectiveness of the linearized complementarity constraint 


signihcant effect as shown in Table In the table, the gap is dehned as 


gap 


(value without constraint) — 
max{|g+p|, 1} 


X 100%. 


5 Conclusion 

In this paper, we have introduced the idea of robust sensitivity analysis for the optimal value 
of LP. In particular, we have discussed the best- and worst-case optimal values under gen¬ 
eral perturbations in the objective coefficients and right-hand sides. We have also presented 
hnite variants that avoid cases of infeasibility and unboundedness. As the involved problems 
are nonconvex and very difficult to solve in general, we have proposed copositive reformu- 
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lations, which provide a theoretical basis for constructing tractable SDP-based relaxations 
that take into account the nature of the uncertainty set, e.g., through RLT and SOC-RLT 
constraints. Numerical experiments have indicated that our approach works very well on 
examples from, and inspired by, the literature. In future research, it would be interesting to 
improve the solution speed of the largest relaxations and to explore the possibility of also 
handling perturbations in the constraint matrix. 
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